Chaos in the Hamiltonian mean field model 
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We study the dynamical properties of the canonical ordered phase of the Hamiltonian mean-field 
(HMF) model, in which N particles, globally-coupled via pairwise attractive interactions, form a 
rotating cluster. Using a combination of numerical and analytical arguments, we first show that 
the largest Lyapunov exponent remains strictly positive in the infinite-size limit, converging to its 
asymptotic value with 1/ln N corrections. We then elucidate the scaling laws ruling the behavior of 
this asymptotic value in the critical region separating the ordered, clustered phase and the disordered 
phase present at high energy densities. We also show that the full spectrum of Lyapunov exponents 
consists of a bulk component converging to the (zero) value taken by a test oscillator forced by the 
mean field, plus subextensive bands of 0(ln N) exponents taking finite values. We finally investigate 
the robustness of these results by studying a "2D" extension of the HMF model where each particle 
is endowed with 4 degrees of freedom, thus allowing the emergence of chaos at the level of single 
particle. Altogether, these results illustrate the subtle effects of global (or long-range) coupling, 
and the importance of the order in which the infinite-time and infinite-size limits are taken: for an 
infinite-size HMF system represented by the Vlasov equation, no chaos is present, while chaos exists 
and subsists for any finite system size. 



PACS numbers: 05.45.-a,05.45.Xt,05.70.Ln,05.90.+m 

I. INTRODUCTION 

In the context of Hamiltonian systems with long-range 
interactions, the Hamiltonian mean-field (HMF) model 
introduced independently in the nineties by Antoni and 
Ruffo 1] and by Kaneko, Konishi and Inagaki 2 became 
the main benchmark for the investigation of thermody- 
namic and dynamical properties of non-additive systems 
[31 2] . The HMF model describes an ensemble of N par- 
ticles moving on a circle, coupled by pairwise (sinusoidal) 
attractive interactions. Each particle can also be seen as 
a pendulum in a fluctuating potential, whose amplitude is 
determined self-consistently and corresponds to the mag- 
netization pp. Detailed studies of the HMF model have 
revealed unusual properties, such as ensemble inequiva- 
lence (associated with the occurrence of negative specific 
heat), long-lived quasi-stationary states, and anomalous 
diffusion [3J 0] . Here, we are interested in the dynami- 
cal properties of the standard (microcanonical) equilib- 
rium phases. Below the critical energy U c — 3/4, the 
HMF system has a finite magnetization (clustered phase) , 
while above U c the magnetization vanishes (homogeneous 
phase). The two regimes are separated by a second or- 
der canonical phase transition. Both in the limit U — > 
and U — y oo the dynamics is integrable: in the former 
case, all particle are trapped in the (harmonic) bottom 
part of the potential well; in the latter, they move freely 
along the circle. At intermediate energies, the (nonlin- 



ear) microcanonical dynamics of a finite system made of 
N particles is characterized by a spectrum of Lyapunov 
exponents (LEs) {A;} with i = 1, . . . , 2N and, due to the 
Hamiltonian structure, A; = — A2w+i-i- 

The thermodynamic limit is, however, a rather intrigu- 
ing subject. For N —> oo, the mean field is constant and 
the evolution of each particle is equivalent to the mo- 
tion of a standard pendulum in a constant gravitational 
field. Accordingly, no chaos but just periodic orbits can 
be generated. This straightforward theoretical prediction 
is consistent with numerical simulations in the homoge- 
neous phase, where it is numerically observed [3J |S] that 
the maximal LE Ai vanishes as A^ 1 / 3 (a result which 
can be easily explained p2 E] by invoking arguments de- 
veloped for products of random matrices [7]). In contrast, 
in the clustered phase, some numerical investigations sug- 
gest that Ai remains finite in the infinite-size limit. These 
findings are consistent with a theoretical study by Firpo 
[8] based on a Riemannian approach, which predicts fi- 
nite Ai values below the transition. However, recently, 
Manos and Ruffo [5] claimed that the A^ -1 / 3 law in the 
homogeneous phase applies also at low energies, specifi- 
cally for U < 0.2, while in the range of 0.2 < U < U c they 
are unable to decide whether the maximal LE vanishes 
or stays finite. Finally, a recent statistical-mechanical 
treatment suggests that the Lyapunov spectrum should 
always converge to zero, but cannot exclude the existence 
of an anomalous subextensive component of strictly pos- 
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itive LEs [TO] . 

In this paper we revisit this issue of the existence and 
nature of chaos in the HMF model. The main part of this 
paper is a study of the largest LE which is split into two 
parts: (i) the analysis of finite-time LEs of a single oscil- 
lator in the fluctuating potential under the assumption 
of a negligible coupling in tangent space; (ii) a careful 
investigation of the effect of the tangent-space coupling. 
The former analysis is justified by the empirical observa- 
tion that the first Lyapunov vector is localized and the 
fact that the influence of a given oscillator on the self- 
consistent mean field, which we call simply the coupling 
strength, decreases as 1/N. We conclude that the single- 
oscillator LE (i.e., the mean of finite-time LEs) cannot 
be larger than 1/lniV, but we also observe that its fluc- 
tuations stay finite in the thermodynamic limit. These 
results are related to the existence of a homoclinic cycle 
connecting the top of the effective potential with itself. 
The following analysis of the coupling in tangent-space 
reveals that even if it is very small it induces a finite 
increase in the LE, which is proportional to the fluc- 
tuations of the single-oscillator LE. This phenomenon, 
that we call "coupling pressure" , is a manifestation of a 
strong sensitivity to coupling which generally arises in 
ensembles of identical weakly- coupled oscillators. It was 
first uncovered in two coupled identical ocillators, where 
it was shown that the maximum LE increases with cou- 
pling strength e by an amount of order 1/| lne| [TT]. The 
same effect was later found in higher dimensional sys- 
tems [T2"HT5] . In the context of globally-coupled systems, 
the coupling pressure is so drastic that it survives even 
though the coupling strength vanishes in the thermody- 
namic limit. We provide a quantitative explanation of 
the effect, by mapping the tangent-space evolution onto 
a stochastic model of sporadically-coupled diffusing par- 
ticles (see Ref. |16j for a preliminary discussion). In the 
HMF model, the effect of the coupling pressure is par- 
ticularly important, since it increases the value of the 
largest LE from zero to a finite number, i.e., it induces 
an instability in a model that would, otherwise, be non- 
chaotic. Altogether, we can summarize our results by 
stating that the infinite-size and the infinite-time limits 
do not commute: taking first the thermodynamic limit, 
the evidence of dynamical instabilities would be lost. An 
indirect confirmation of the theoretical approach comes 
from the localization of the Lyapunov vector, that is con- 
firmed by our numerical simulations. 

Our investigation of the largest LE in the ordered 
phase of the simple HMF model is completed by the study 
of a number of related points: first, we numerically study 
the largest LE in the vicinity of the critical energy value 
U c . We find that Ai goes to zero for U — > U c from below 
and account for the observed scaling behavior. 

Next we address the problem of the shape of the en- 
tire Lyapunov spectrum. We find that several exponents 
(in addition to the largest) stay positive, but their num- 
ber is non-extensive (i.e., it grows slower than linearly 
with N). The results are consistent with the theory in 



Ref. [TO] where it was predicted that "with measure one" 
the spectrum should be equal to zero. 

Finally, we discuss a 2D generalization of the HMF 
model jTTl |T5] , to test the general validity of our theoreti- 
cal and numerical findings. In the 2D model, each oscilla- 
tor is composed of four variables and thus can be chaotic 
without taking the coupling pressure into account. We 
find nevertheless the same size-dependence of the largest 
LE as in the standard HMF model. We also investi- 
gate the full Lyapunov spectrum in this case and argue 
to what extent the argument developed for the standard 
HMF can be extended here. 

This paper is organized as follows: The HMF model 
is introduced in Sec. |Hj together with a careful discus- 
sion of its equilibrium properties. This is necessary to 
collect proper information on finite-size effects that is 
crucial for a correct development of our theoretical argu- 
ments. Section IIHI is devoted to a critical discussion of 
the numerical results for different system sizes and differ- 
ent energy values. There, we illustrate some of the issues 
that hindered the interpretation of the numerical results. 
Section IIVI is devoted to a detailed characterization of 
the evolution in the tangent space of a single particle in 
a self-consistent mean field. In particular, we introduce a 
finite-time LE and discuss its dependence on the energy 
and the number of particles. The effect of the coupling 
is discussed in Sec. [V] where we first introduce a simpli- 
fied model and test the correctness of our solution. The 
application to the HMF model is analyzed in the second 
part of the section. The scaling of the largest LE at the 
transition energy U c is derived in Sec. |VI|and compared 
with numerics. The structure of the Lyapunov spectra is 



analyzed in Sec. VII In Sec. VIII we deal with the 2D 
generalization of the HMF model. A brief summary of 
the results and a discussion of open questions are finally 
reported in Sec. llX) 



II. 



THE HAMILTONIAN MEAN FIELD 
MODEL: EQUILIBRIUM RESULTS 



The HMF model was derived from a one-dimensional 
self-gravitating model, by truncating the Fourier expan- 
sion of the gravitational potential to its first term [1] . It 
consists of N unit-mass particles that move on a circle 
under their mutual attraction. The dynamics of the N 
particles is ruled by the Hamiltonian |19j 
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where Oi and pi denote particle positions (angles) and 
velocities. The resulting equations of motion write 



h) = Msin( 
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where M is the magnetization and 
phase, defined by 



N ^ 



the associated it is equal to half of it, namely 



(3) 



M measures the degree of clusterization and plays the 
role of an order parameter |20j . Depending on the en- 
ergy density U = H/N, the system can show two differ- 
ent thermodynamic phases, separated by a second-order 
transition: (i) the clustered, ordered phase, characterized 
by a finite magnetization (for U < U c = 3/4) ; (ii) the 
homogeneous phase, characterized by a vanishing magne- 
tization (for U > U c ). In the following, we limit ourselves 
to the clustered phase U < U c (T < T c ). 

All the reported simulations have been performed 
within the microcanonical ensemble by implementing 
symplectic integration schemes, typically a 4-th order 
McLahlan-Atela algorithm [21] with integration time step 
dt = 0.05 or 0.1. This choice ensures an energy conser- 
vation with a relative precision of the order of 10~ 10 to 
10" u . 

Initial phases and momenta have been typically drawn 
from the invariant equilibrium distribution discussed in 
the next subsection. We have also occasionally compared 
the results with those obtained for different choices of 
initial conditions, namely: (i) zero phases and a Gaus- 
sian distribution of the momenta; (ii) a uniform distri- 
bution of the phases and a Gaussian distribution of the 
momenta. A transient (typically 5x 10 3 7V time units) has 
been discarded, before starting the computation of any 
equilibrium quantity. Finally the typical transient time 
for the evolution in tangent space lies between 4 x 10 5 
and 4 x 10 6 time units, while the typical integration time 
lies between 4 x 10 6 and 10 7 time units. 



A. Equilibrium distribution of single oscillator 
energy 

From the point of view of a single oscillator, the evolu- 
tion equation in Eq. ([2]), at finite system sizes, is equiv- 
alent to that of a pendulum in a noisy environment, the 
noise being the result of the statistical fluctuations of the 
magnetization. In the thermodynamic limit, M and 4> 
are strictly constant and, as a result, the single oscillator 
energy 
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is strictly conserved. Notice that, in Eq. Q, the (arbi- 
trary) zero level of the potential energy is shifted by 
I'm = M — 1 in order that the ground-state energy of the 
single oscillator is always zero for any value of the magne- 
tization M. This is the convention adopted throughout 
the paper Note also that the total potential energy 
is not given by the sum of the single potential terms, but 
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the reason being that each term would otherwise be 
counted twice. 

The equilibrium distribution of the single-oscillator en- 
ergies is 
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Q(P,0,T), (6) 



where Q{p 1 9,T) is the Gibbs-Boltzmann distribution 
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Q(p,8,T) = Cexp 
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with a suitable normalization constant C, the unit Boltz- 
mann constant, and the temperature T given by [TJ [31 2] 



U= T ~ + \{1-M>). 



As a result, one finds that 
P{h,T)-- ( 



-h/T 



2MJo y/y(h/M-y)(2-y) 



(8) 



dy, (9) 



where yo — h/M if h/M < 2 and yo = 2 otherwise. The 
integrand has two (integrable) square-root singularities 
at both ends of the integration interval for all energy 
values, except for h — 2M, in which case the singular- 
ity is hyperbolic, and this indicates a logarithmic diver- 
gence of the integral. Note that the divergence lies at the 
separatrix e s of the single-particle effective Hamiltonian 
Q. This equilibrium energy distribution ([9]) is plotted in 
Fig. [lja) for three different energy densities (U — 0.15, 
0.2 and 0.5). 



B. Magnetization 

In this section we discuss the behavior of the magne- 
tization in the clustered phase and at the critical energy 
U c for finite N. For large but finite N, the magnetiza- 
tion is affected by statistical fluctuations and, as a re- 
sult, the oscillator energies diffuse, albeit very slowly. 
Simple statistical arguments suggest that the absolute 
value of the magnetization M{t) fluctuates around its 
mean field value with an amplitude that should scale 
as 1/yN. We numerically checked this conjecture by 
measuring the power spectrum S(f) of the magnetiza- 
tion M(t) for different system sizes. In Fig. [ljb), we see 
that most of the power is concentrated in a broad peak 
around / = 0.16 (for U = 0.5) and the peak power scales 
as I/AT, in agreement with the l/y/~N amplitude of the 
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FIG. 1: (color online), (a) Equilibrium energy distribution 
[Eq. Q] for three different values of the internal energy: U = 
0.1, U = 0.2, and U = 0.5 (from bottom to top). The value of 
T in Eq. (|9j) is computed by Eq. Q . Recall the shifted energy 
axis, Eq. d4||, for h. (b) Power spectrum S(f), multiplied 
by the size N, for U = 0.5 and three different system sizes, 
N = 250 (solid black line), 500 (dotted dashed red line), and 
1000 (dashed green line). 



fluctuations. We also verified that the power contained 
in the low-frequency peak increases upon increasing the 
system size. The peak location should be related to some 
characteristic timescale of the dynamics, though we have 
no precise hints about its origin. 

In contrast, the motion of the global phase is deter- 
mined by the oscillators with single-particle energies h 
larger than the separatrix energy e s pQ. They continue 
rotating either clockwise or counterclockwise according 
to their momentum p. The slow energy diffusion of in- 
dividual oscillators implies that those with low energies 
wander and can reach h > e s , "randomly" picking the 
rotation direction, and stay in this high-energy state for 
some time until they eventually go back to the "bounded" 
state with h < e s . The numbers of particles rotating 
clockwise and counterclockwise are on average equal to 
one another, but, because of statistical fluctuations, the 
instantaneous fractions of the populations typically dif- 
fer by a quantity of order 1 /y/~N. Because of momentum 
conservation, the phase <f>^ [23 of the global magnetiza- 
tion exhibits a net drift with an average angular velocity 
w ~ l/y/~N PQ, which has also been verified numerically 
(not shown). Over long time scales, the sign of the ve- 
locity changes since the fluctuations will invert the pre- 
dominance of clockwise/counterclockwise rotating parti- 
cles. As a consequence, we expect the global phase to 
exhibit a crossover from drifting to diffusive motion over 
a crossover timescale T<jiff. Numerical simulations (see 
Fig. [2]) clearly show that the crossover time diverges in 
the thermodynamic limit as idiff ~ N. 

Finally, we discuss the equilibrium behavior in prox- 
imity of the critical energy. In the thermodynamic limit, 
the magnetization M obeys the usual mean-field behav- 
ior, i.e., M~\U- £4 1 1/2 for U < U c J2i]. Determining 
M for finite N is, however, a delicate problem which re- 
quires a careful treatment of finite-size fluctuations [2~11 - 
121)] . The correct solution can be found by taking into 



FIG. 2: (color online), (a) Mean-square displacement of the 
global magnetization phase A<^) 2 = ((4>(t) — (f)(0)) 2 } vs. time. 
Mean-square displacements have been computed by averag- 
ing over time span of order 1 x 10 7 with sliding windows 
of duration 160000 — 320000 for three different system sizes: 
namely, N = 1000 (black dots), N = 4000 (red squares) and 
N — 10000 (green diamonds), (b) Same quantities as in (a) 
but rescaled to better highlight the crossover from ballistic 
to diffusive behavior. Time has been rescaled by system size, 
t' = t/N, and mean-square displacements by the rescaled 
time, A(f) 2 — > A(f) 2 /t' . The dashed lines mark the linear bal- 
listic growth, and the vertical dot dashed (blue) line highlights 
the beginning of deviations from linear growth at the rescaled 
crossover time Tam/N. 



account the law of large numbers in the self-consistency 
equation of the mean field argument. The magnetization 
can be expressed within the canonical ensemble formula- 
tion as 



M + SM = 



/oo p2tt 
dp / d9e* e c- n / T 
-00 Jo 



(10) 



where H — p 2 /2 + 1 — Mcos# (here the absolute scale 
is chosen for the energy), Z = dp J Q * d6e~ n / T , 
and I • I denotes the modulus of the complex number. 
The second term in the l.h.s. represents an unavoidable 
finite- N correction that we assume to be in the order of 
SM ~ 0(M/y/~N). A straightforward calculation (see 
also Ref. [1]) leads to the self-consistency equation 
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where I n (z) is the first-kind modified Bcssel function of 
order n. Near the critical point, Eq. ( 11 ) can be expanded 
for small M, yielding 
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In the infinite size limit, SM — 
expected mean- field result 0], 
T < T c = 1/2. For finite sizes N, Eq. |12j) predicts the 
following scaling, 



0, and this gives the 
M - (T c - T) 1 / 2 for 



M(U,N) - N~P /u F 



((T c - T)N 



for T<T C , 
(13) 
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FIG. 3: (color online). Critical scaling of the magnetization 
M. (a) Size dependence of the magnetization at the critical 
point, M(U C ,N). (b) Magnetization M(U,N) as a function 
of the energy difference AU = U - U c for N = 100, 400, 1000, 
and 4000 (from top to bottom), (c) Same data as the panel 
(b) with rescaled axes. 



with f3 — 1/2, v — 2, and a scaling function F(z). Using 
the energy-temperature relation ^ [27] . we can rewrite 
Eq. (pi) as 



M{U, N) ~ N~ pfv G ((£/ c - U)N 1/v ^ 



for U < U c , 
(14) 

with another scaling function G(z). This expression 
accounts for the critical decay of the magnetization 
M(U C ,N) ~ N- 1 / 4 found in Fig. [^a). Equation §H§ 
can be further checked by rescaling the magnetization 
M(U.N) off criticality; plotting MN@/ V against (U - 
U^N 1 /" , we confirm that the data shown in Fig.[3jb) col- 
lapse reasonably well onto a single curve, G(z) [Fig.pjjc)]. 
It is worth noticing that the observed finite-size scaling 
M(U C ,N) ~ iV -1 / 4 for the magnetization was reported 
for the first time for the mean field version of the Ising 
and Heisenberg model in Ref. [35] . The obtained value of 
the critical exponent v — 2 is the one found for dissipative 
noisy phase oscillators [25], but is different from that of 
the Kuramoto model, i.e. deterministic phase oscillators 
with random frequencies, v = 5/2 |26) . The value v = 2 
is also the one expected from a simple dimensional analy- 
sis, v — d c , where %p is the usual correlation- length 
exponent in the mean-held limit and d c is the upper crit- 
ical dimension |29j. In our case, d c — 4 and ^mf = 1/2, 
which further confirms the analogy of the HMF model 
with the mean held XY Heisenberg model [T|. 



III. LYAPUNOV CHARACTERIZATION OF 
THE DYNAMICS 

In order to characterize the dynamics of the system, 
we estimate the LEs by following the dynamical evolu- 
tion in tangent space of a vector v — {69i, 5pi}i = i....,N of 
infinitesimal perturbations, 



59t = Sp t , 

Spi = —M cos(0 - ed&Oi + ■= ^2 ' 



(15) 




FIG. 4: (color online). Largest LE Ai of the HMF model, 
(a) Ai versus the system size iV for U = 0.5 (black cir- 
cles) and U — 0.7 (red squares). The dashed lines show the 
quadratic extrapolations to the asymptotic values. (b,c) Time 
series of the finite time exponent A r for N — 100, 400, 1000 
with U — 0.25 (panel (b) from top to bottom) and for 
U = 0.20,0.25,0.30 with N = 400 (panel (c) from top to 
bottom). The laminar and burst states have different values 
of the time-averaged LE as shown in Fig. |6|c) below, (d) 
Time fraction in the burst state Fb urst as a function of the 
size iV for U = 0.25 (black circles), 0.30 (red squares), and 
0.35 (green diamonds). 



and by orthonormalizing the resulting vectors at proper 
times [3D] • In order to probe large system sizes up to N = 
10 6 , we have implemented a highly parallelized version of 
the Gram-Schmidt algorithm. 

Most of the numerical simulations have been performed 
for U — 0.7 (which corresponds to a magnetization 
M ss 0.281 and a temperature T 0.479, measured at 
N > 10 5 ) and U = 0.5 (M 0.621 and T 0.386). 
Both parameter values are sufficiently away from the crit- 
ical point U c as well as from the zero-temperature limit; 
otherwise the asymptotic behavior of the Lyapunov ex- 
ponent would be masked by severe finite-size effects (see 
below) . 

According to a theoretical argument briefly sketched 
in Ref. [TH], where we showed that, in globally-coupled 
dissipative systems, the leading finite-size corrections to 
the asymptotic value of the largest LE are polynomial in 
1/lnTV, we find it convenient to investigate the finite-size 
dependence of Ai by plotting it as a function of 1/lniV. 
The data reported in Fig. Qa) are indeed consistent with 
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a logarithmic dependence 
Ai = A c 
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In 2 N 
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for large N, especially for U — 0.7 (red squares). De- 
viations from the 1/ In iV behavior are visible for small 
N, but taking a quadratic correction into account is 
sufficient to describe perfectly all system sizes studied 
(Fig. Qa), red dashed line). This quadratic correction 
is stronger for U = 0.5, which is an incipient evidence 
of the convergence problems that arise at small energies 
(see below). Because of this slow but significant size- 
dependence, the extrapolated values of the maximum LE 
in the thermodynamic limit are smaller than the typical 
values reported in the literature (see, e.g., Ref. [3]), but 
are neverthless clearly different from zero: Aco = 0.056(6) 
atU = 0.5 and A^ = 0.046(3) at U = 0.7, both obtained 
by using the quadratic ansatz described above. 

As briefly mentioned in the introduction, a number 
of numerical studies have reported contradicting conclu- 
sions about the largest LE in the clustered phase: while 
most of them claimed, qualitatively, no or weak size- 
dependence and thus strictly positive asymptotic values 
of the maximal LE [3J [5] , Manos and Ruffo [3] reported 
a power law decay A ~ TV -1 / 3 , although their simula- 
tions were performed for substantially lower energy den- 
sities (U = 0.1). Our own numerical simulations (not 
shown) performed at U = 0.1 indeed confirm the power- 
law decay at least up to N = 10 6 . In order to shed 
some light on the possible existence of two qualitatively 
different phases, we scanned intermediate energy levels 
in the interval U £ (0.1,0.5). These simulations (see 
below) revealed the presence of strong intermittent be- 
havior, which make practically impossible to determine 
a reliable value of the LE, in particular for large system 
sizes. The phenomenon is better illustrated by studying 
the finite time LE 



t \\v(t-r) 



(17) 



fixing r = 2 [31]. Fig. Qb,c) reveals irregular jumps of 
A T (<) between two clearly different states: (i) a laminar 
one, where X T (t) stays near zero; (ii) bursts, where A T (i) 
fluctuates much more strongly. Upon comparing simula- 
tions performed for different sizes [Fig.j4j[b)] and different 
energy densities [Fig. pffc)], we see that the frequency of 
the bursts grows both with N and U. This is quantita- 
tively shown in Fig. Qd) . 

In order to understand the origin of this intermittent 
behavior, we analyze the structure of the (first) Lyapunov 
vector. The Lyapunov vector is the quantity associated 
with each LE and indicates the direction of the infinitesi- 
mal perturbations growing at the rate of the correspond- 
ing LE. It is defined as a function of the phase-space 
point and turns out to be a useful tool to characterize 
statistical properties of large dynamical systems [3"2"H3"4"] . 
Here, it is convenient to introduce the squared amplitude 




FIG. 5: (color online). Localization of the Lyapunov vector, 
(a) a, vs i for N = 4000 (black solid line), N = 10000 (red 
dotted dashed line), and N = 20000 (green dashed line) at 
U = 0.5. The dotted black line marks a decay as 1/i. (b) 
Average amplitude of the vector components as a function of 
the corresponding single-oscillator energy, a(h), in a system 
of N = 10 oscillators at U — 0.5. The vertical dashed red 
line marks the single oscillator separatrix energy e s = 2M. 



of the normalized vector components for each oscillator, 

A 4 = 8% + 8tf, (18) 

and to consider its time average (Ai) (here and in the fol- 
lowing, angular brackets denote time averages), to have 
a statistically reliable quantity. 

In homogeneous globally-coupled systems, any order- 
ing of the oscillators is equally meaningful, as they are 
equivalent to one another. In Fig. [5ja) , we plot the am- 
plitude ai — \J (Ai) versus its rank (i.e., we arrange 
the oscillators according to in decreasing order) for 
N = 4000, 10000, and 20000. Our data show that the 
Lyapunov vector is approximately localized as 1/i (as 
indicated by the dotted black line), that is, the pertur- 
bation is concentrated in a few components. In Fig.[HJb), 
the data is organized in a different way. The oscillators 
are grouped according to their energy and the pertur- 
bation amplitude is averaged over all oscillators in the 
interval [h, h + dh], to obtain a(h). The results reported 
in Fig. [5jb) indicate that the vector component is sub- 
stantially larger when the energy of the corresponding 
oscillator is close to that of the separatrix. 

In order to further clarify the relationship between lo- 
calization and energy, it is instructive to monitor the in- 
stantaneous degree of localization of the Lyapunov vec- 
tor, by estimating the inverse participation ratio |35j 



Yi 



(19) 



By construction, 1/N < Y 2 < 1. The larger is Y 2 , the 
more the vector is localized; Y 2 = 1 denotes complete lo- 
calization on a single oscillator, while Y 2 = 1/N indicates 
a completely delocalized vector with equal components. 

In Fig. [6]ja), we plot the time evolution of the finite- 
time LE A r , of the inverse participation ratio Y 2 , and of 
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FIG. 6: (color online). Intermittency for U = 0.25 and 
N = 400. (a) Time series of the finite-time LE A T (time 
window t — 2), the instantaneous value Y2 of the inverse par- 
ticipation ratio, and the energy difference hyi — e s between 
the dominating oscillator (the one with the largest amplitude 
Ai of the Lyapunov vector component) and the separatrix. 
(b) Density plot with respect to /im — e a and Yi. The color 
code shows the frequency in the logarithmic scale. The black 
indicates null density. Two peaks corresponding to the burst 
and the non-burst states are clearly visible, (c) Lyapunov 
exponent Ai (black circles) versus N for U = 0.25. The con- 
ditioned Lyapunov exponents (see text) are also shown for the 
burst state (red triangles) and for the non-burst state (green 
diamonds) . 



the energy hja of the oscillator with the largest amplitude 
Ai in the Lyapunov vector components. The data refer to 
a small system (TV = 400) with energy density U = 0.25. 
The three temporal traces reveal a strong correlation be- 
tween the occurrence of the bursts in the finite-time LE, 
a stronger localization and the closeness of the energy to 
that of the separatrix. A more quantitative characteriza- 
tion of the connection between the inverse participation 
ratio and the energy is presented in Fig. [6jb) , where the 
color code indicates the probability to observe a given 
pair of values (Y 2 , Km — e s ). Altogether, the data plotted 
in Fig. J6ja,b) confirm the intermittency between the two 
distinct states: (i) the laminar state is characterized by a 
less fluctuating finite-time LE, a weak localization of the 
Lyapunov vector, and single-oscillator energies far from 
the separatrix; (ii) the burst state by large fluctuations 
of the finite-time LE and a strong localization of the Lya- 
punov vector around an oscillator lying very close to the 
separatrix. 

In Fig. |6][c) we compare the value of the true, time- 
averaged LE (circles) with the averages restricted to the 
bursts (triangles) and the laminar state (diamonds) for 



U = 0.25. We see that upon increasing the system size, 
the "burst" LE tends to converge towards the true LE. 
This reflects the fact that the laminar state tends to dis- 
appear for N — > 00 and the "laminar" LE remains quite 
small. 

Therefore, we conclude this numerical analysis by 
noticing that the observed value of the LE depends 
strongly on whether there is at least one oscillator whose 
energy is sufficiently close to the energy e s of the sep- 
aratrix. If, for any reason, no oscillator has an energy 
hi close enough to e s , the laminar contribution domi- 
nates. Altogether, our analysis suggests that a truly 
asymptotic behavior is observed only if (on average) 
at least one oscillator has an energy sufficiently close 
to that of the separatrix. The minimal number N m 
ensuring this condition can be estimated by imposing 
N m P(e s ,T)Sh = 1 with a suitable width Sh of the en- 
ergy window. It grows quickly with decreasing U, since, 
for low energies, the energy distribution is approximately 
exponential, P(h,T) rj exp(-h/T) with T m 2U [see 
Eq. (J8J], while the separatrix energy is practically con- 
stant, e s = 2 AI m 2. By referring to the theoretical 
expression in Eq. ^ and assuming %/iV (see the 

next section for a justification), we find that N m m 10 13 , 
10 5 , and 10 for U = 0.1, 0.2, and 0.5, respectively (these 
are the energy values considered in Fig. [I]). It is clear 
that for U — 0.1 there is no hope to reach the asymptotic 
regime in numerical simulations with the currently avail- 
able machines and, in particular, that the Ai ~ TV -1 / 3 
scaling found by Manos and Ruffb [5j characterizes only 
the laminar state, which is not the asymptotic state of 
the system. 

Although our numerical results strongly support the 
strictly positive asymptotic value of the maximal LE, it 
is not clear how this behavior is connected to the pres- 
ence of oscillators in the vicinity of the separatrix. The 
following two sections are devoted to clarifying this point. 



IV. SINGLE OSCILLATOR ANALYSIS 

In this section we analyze the behavior of the Lya- 
punov exponent and vector, by neglecting the coupling 
term in the tangent space [i.e., the sum in Eq. ( 15 1]. This 
assumption is tantamount to studying a single oscillator 
forced by the field Mjv(i)e^ JV (*- ) , which is generated self- 
consistently by an ensemble of N globally-coupled oscilla- 
tors. The evolution is ruled by the effective Hamiltonian 
Q and is thereby described by the equation 



= p 

p = Mpf sm(4>N — ' 
which, in the tangent space, becomes 
58 = 5p 

Sp = —Mpj cos(4>n — 



(20) 



(21) 



with no contribution from the coupling with the other 
oscillators. As already noted, the two observables Mm 
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FIG. 7: (color online). Single-oscillator LE versus the energy 
of the single particle, forced by TV globally-coupled oscillators 
at U = 0.5. In panel (a) the solid curves correspond to TV = 
100, 1000, 4000, and 10000 (from top to bottom); the three 
rescaled curves in panel (b) correspond to TV = 1000, 4000, 
and 10000. The curves overlap except in the peak region, or 
near the separatrix energy. The vertical dashed lines indicate 
the energy e s of the separatrix. Panel (c) shows a close-up of 
the peak at h = e s , the LE being scaled by the maximum value 
Am of each curve and plotted against rescaled single particle 
energies (crosses, triangles, diamonds, and circles correspond 
to TV = 100, 1000, 4000, and 10000, respectively). 



and 4>n are strictly constant in the thermodynamic limit 
TV — » oo. This means that the corresponding LE is ex- 
pected to be equal to that of a standard pendulum, i.e., 
zero. 

In the following, we investigate the size-dependence of 
the maximum LE of the single forced oscillator, by in- 
troducing an energy-dependent single-particle LE Xo(h). 
Since a meaningful definition of a LE involves the infinite- 
time limit, while energy is conserved only during finite 
times at finite TV, we introduce sporadic small correc- 
tions to prevent the trajectory from diffusing away from 
the prefixed energy shell h. This is achieved by rescaling 
the kinetic energy, or the particle velocity, each time the 
trajectory passes through the point of minimal potential 
energy, without adjusting the potential energy. 

The numerical results reported in Fig. [7] confirm that 
the LE decreases with increasing TV as expected, but it 
also displays a strong dependence on the energy. The 
peak is centered at the energy e s of the separatrix. In 
panel (b) we see that everywhere except in the peak area, 
Ao scales as TV^ 1 / 3 . This behavior can be understood by 
invoking known results for random symplectic matrices, 
similarly to the maximal LE of the full system in the ho- 
mogeneous phase. It is known [7] that approximating the 
tangent-space dynamics with a product of independent 
symplectic random matrices with zero-mean disorder of 
amplitude r\ gives rise to a positive Lyapunov exponent 
which scales as r; 2 / 3 . In our setup, the statistical fluctu- 
ations of the collective magnetization scale as 1/ \/TV and 
play the role of the disorder. As a result, r\ ~ 1/yN and 
this explains the —1/3 scaling clearly seen in Fig. [7jb). 

The same argument, however, does not apply to the os- 




FIG. 8: (color online) . Comparison between the maximum LE 
Am of a single forced oscillator (black squares) at separatrix 
energy e s and the first LE Ai of the full HMF model for U = 
0.5 (red circles). The dashed lines highlight the quadratic 
behavior Ai, A M ~ 60 + bi/ln/V + b 2 /(hiTV) 2 with b = for 
Am- (a) The LEs versus 1/lnTV. (b) The LEs are multiplied 
by In TV. The absence of the constant term 60 is confirmed by 
the linearly arranged symbols for Am- 



cillator with energy near e s , because here the instability 
is rather due to the separatrix. This results in a peak in 
\o(h) at the separatrix energy, which does not decay as 
TV -1 / 3 . We investigate the scaling behavior of its width, 
by plotting in Fig. [7^c) the LE normalized by its max- 
imum value Am (for any given TV) with a rescaled axis 
(h — e^N 1 / 2 . The nice overlap of the curves obtained 
at different system sizes shows that the width decreases 
as l/\/TV. This indicates that the anomalous behavior is 
exhibited by 0{\N) oscillators located within the range 
of the separatrix energy. These are consequences of the 
0(1/ y/N) fluctuations of the magnetization. 

Finally, in Fig. [8] we plotted Am for different values of 
TV (see black squares). The data shows that Am scales 
as 1/lnTV for large TV. This behavior can be understood 
by introducing a suitable symbolic dynamics. The main 
source of uncertainty (and thus of entropy) is associated 
to the binary "choice" made by the oscillator on reaching 
the top of the potential, between the option to return to 
the same side or to pass it. Accordingly, we expect the 
metric entropy to be K In 2 jt s , where t s is the return 
time to the saddle [36] . The return time can be estimated 
as the time needed to amplify a distance from the sad- 
dle, in the order of the noise amplitude 1/vTV, to a value 
of order 1. Since the separation rate from the saddle is 
finite in the thermodynamic limit, the condition reads 
exp(t s )/y/~N 1. Accordingly, t s ~ In TV and therefore 
K « 1/lnTV. Since the metric entropy is generically es- 
timated by the sum of the positive LEs, which is simply 
equal to the sole positive LE in our case, we can finally 
conclude that Am ~ 1/lnTV, too. This prediction is con- 
firmed in Fig. [8] with a quadratic correction for finite TV, 
i.e., A M ~ 61/ In TV + 6 2 /(lnTV) 2 + . . . . Note that, in con- 
trast to the first LE Ai of the full system (red circles in 
Fig.[8|, Am vanishes in the infinite-size limit as it should 
be (black squares). 
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By further comparing the single-oscillator LE Am with 
the first LE Ai of the full system (see red circles in Fig. [8]), 
we see that such a single oscillator contribution shows a 
size-dependence similar to that of the full-system LE. 
However, for increasing N there is no evidence that the 
gap is going to close. Indeed, Fig.[8jb) shows that Ai In A 
diverges for A — > oo, suggesting that in the thermody- 
namic limit the relevant, non vanishing contribution to 
the first LE Ai arises from the coupling terms. It is there- 
fore necessary to consider more carefully the whole evo- 
lution equation in tangent space and the role played in 
this context by the fluctuations of the finite-time LE. 
A simplified argument can be put forward to infer the 
asymptotic behavior of the fluctuations of the finite-time 
LE for a forced oscillator near the separatrix energy. 
Whenever the oscillator passes near the saddle, its growth 
rate is always of order 1, positive or negative, irrespec- 
tive of the system size. This implies that the amplitude of 
the finite-time LE fluctuations near the separatrix energy 
should remain positive in the thermodynamic limit, pos- 
sibly with logarithmic finite-size corrections. In the next 
section, we will analyze these fluctuations more closely 
in order to build a simplified model for the fully coupled 
tangent space dynamics which will allow us to establish 
a finite lower bound for the largest LE. 



V. A SIMPLIFIED MODEL FOR THE 
COUPLING PRESSURE 

In the context of globally-coupled dissipative systems, 
we recently showed that global coupling may induce an 
increase of the first LE with respect to the single-unit 
exponent |16) . Here, we refine such argument for the 
HMF context. 

The single oscillator approximation discussed in the 
previous section consists in disregarding the contribu- 
tion of the coupling term appearing in the second line 
of Eq. (15). In fact, as we explain below, this is what 



happens to the tangent-space evolution of the full sys- 
tem for most of the time, because of the localization 
of the Lyapunov vector (see Fig. [5]) and the 1/N nor- 
malization in front of the coupling term. The Lyapunov 
vector components then evolve independently. In partic- 
ular, the logarithms of the amplitudes (xi = In y/A~i) can 
be regarded as Brownian particles with a drift velocity 
given by the single-particle LE and an effective diffusion 
constant that measures the fluctuations of the LE itself. 
The analysis carried out in the previous section suggests 
that the oscillators should be classified into two groups: 
(i) a small fraction which lies close to the separatrix and 
is characterized by a LE of order 1/ In A and non zero 
fluctuations of the finite time LE; (ii) the vast major- 
ity, characterized by a LE of order A -1 / 3 and vanishing 
fluctuations. The 0{l/y/N) fluctuations of the magne- 
tization then suggest that the population ratio between 
the particles close to the separatrix and the remaining 
population vanishes in the thermodynamic limit, possi- 
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FIG. 9: (color online). Schematic representation of the toy 
model for the coupling pressure (see text). Point-like particles 
sit at the logarithmic coordinates Xi — In y/Ai and are either 
diffusive (red circles) or non diffusive (blue). Diffusive parti- 
cles turn non diffusive with a rate ai and non diffusive ones 
start to diffuse with a rate ~ otx/yN, as the correspond- 
ing oscillators approach and leave the separatrix, respectively. 
The coupling is zero as long as particles are no farther than 
In N from the rightmost particle, but otherwise acts as a bar- 
rier, preventing particles from being farther than In TV from 
the rightmost one. The net effect is a drift to the right. 



bly as 1/vA. Moreover, energy diffusion induces (slow) 
exchanges between the two families. 
We now discuss how the coupling modifies the single os- 
cillator evolution. The localization of the Lyapunov vec- 
tor indicates that the coupling term (the sum in the r.h.s. 
of Eq. ( [15] )) is of the order of S9 m /N, where m labels 
the oscillator where the vector is localized. Therefore, 
because of the 1/N factor, the mth oscillator component 
only weakly affects the oscillator at stake, 88i, and so does 
the coupling term. However, the opposite is true when 
\S9i\ <C \69 m \/N (notice that this is possible, since the 
various components evolve independently of each other, 
and their logarithms diffuse away). In this latter case, the 
evolution is dominated by the coupling and the net result 
is that such extremely small components become of the 
order of the coupling term. In terms of the logarithmic 
coordinates Xi, the effect of the coupling can be schema- 
tized by a barrier sitting at a; m ; n = x max — In N (where 
i maK labels the rightmost particle, that is the largest vec- 
tor coponent) which prevents any interparticle distance 
from being larger than In N. 

Altogether, we propose a simplified model of two pop- 
ulations of "particles" , as sketched in Fig. [9j The par- 
ticles in the first group (red in Fig. [9| show the biased 
Brownian motion, while those in the other group (blue) 
stay quiescent. Each particle evolves independently of 
the others until it lies at a distance larger than In A to 
the left of the rightmost particle (where the Lyapunov 
vector is localized), in which case it is instantaneously 
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pushed forward by the coupling to restore the maximal 
allowed distance (drawn by the gray box of size In N in 
Fig. [9]) . A precise formalization of the model requires the 
following additional ingredients: drift velocities and dif- 
fusion coefficients of the two populations and the mutual 
transition rates. For the drift, we assume that both pop- 
ulations are characterized by a zero velocity (zero LE). 
Since our goal is to explain the origin of a strictly positive 
LE in the thermodynamic limit, we believe that neither a 

1 /TV 1 / 3 nor a 1/ In N LE can eventually provide a leading 
contribution and thereby set both to be zero. As for the 
diffusion coefficient of the Brownian particles, we assume 
a finite value D s for the first population (that correspond- 
ing to the finite-time LE fluctuations of the oscillators in 
the vicinity of the separatrix). In contrast, we assume a 
zero diffusion coefficient for the second population, as it 
is negligible for large sizes. As for the transition rates oi 
and a.2 from the diffusing to the still population and vice 
versa, respectively, on the basis of the numerical obser- 
vation and the theoretical argument on the ratio of the 
two populations, we assume that the ratio 0:2/0:1 van- 
ishes in the thermodynamic limit as \j\J~N . Thus, we 
are left with three independent parameters: the diffusion 
coefficient D Sl the transition rate from the diffusing to 
the still population ai, plus a small parameter 02 of the 
order of ai/yN. 

It is convenient to introduce the probability density 
Pj(x,t) for a particle of the jth population (j = 1 and 

2 referring to the diffusing and still populations, respec- 
tively) to be at position x. If both populations move 
with a positive velocity w, a positive LE spontaneously 
emerges in the system of interacting particles. It is con- 
venient to study the problem in a frame moving with a 
velocity v, since then one has to look for a stationary 
solution. In this frame, the evolution equation reads, 



Ik 

dP 2 
Ft 



v di\ + p 

dx 2 dx 2 1 1 
dP 

v—— + u x P x - a 2 P 2 , 
ox 



a 2 P 2 



where x £ [0, 00], with a reflecting boundary at x = 0. 
Equation (22) describes an ensemble of stochastic parti- 



cles that move along a tilted plane, which corresponds 
to the velocity v of the comoving frame, and have two 
possible internal states, one characterized by a finite dif- 
fusion D s and the second one by a zero diffusion, so that 
the particles in the second state simply move toward the 
reflecting barrier at x — 0. The diffusive dynamics com- 
petes with the time scales set by the transition rates be- 
tween the two populations. In particular, if D s is finite 
and oi and 02 are sufficiently small, particles in the non 
diffusing populations will tend to accumulate in x = 0. 
Therefore, we look for a general stationary solution of 
Eq. (1221) of the form 



-JX 



c S(x) 



c 2 e 



-■yx 



(23) 



with some 7 > and Dirac's delta S(x). The two proba- 
bility densities must obey particle conservation 



(Pi(aO + P 2 (a:)) dx = l, 



and the population equilibrium condition 



01 / Pi(x) dx 
'0 



02 / P 2 (x)dx. 
10 



(24) 



(25) 



Substituting the Ansatz (23) into Eq. (22), we obtain the 
stationary conditions in the bulk (x 7^ 0), 



= — WC17 + 
= — vc 2 j + 
which yields (for d/0) 



D. 



s 2 

7 ci 



2 

OiCi 



a 2 c 2 , 



o 2 c 2 , 



C2 



;d s1 



(26) 



(27) 



and 



2 7 v 2 + [2(ai + a 2 ) - D sl 2 } v - D s a 2 -f = 0. (28) 

This can be solved for v, choosing the physically mean- 
ingful positive solution. By recalling that 02 ~ oi/V^/V, 
we can expand in terms of o 2 , 




+ O (o 2 ) 



2ai-_D a7 



O(at) 



if (o! + o 2 ) < Z) s 7 2 /2, 
if (oi + o 2 ) > D sl 2 /2. 

(29) 

Note that the velocity (i.e., the LE) is strictly positive 
when diffusion dominates over the interstate transitions. 



From Eqs. (|23}-([25} it also follows 
o 2 /oi 



ci = 7 



1 + 02/01 



o 2 
' — . 

01 ' 



(30) 



and together with Eqs. (27) and (29), we find that the 



(22) coefficient of Dirac's delta has a finite amplitude, 



(jo 
2 



l-O (a 2 ) if (o!+o 2 ) <D sl 2 /2, 
%£-<D(a 2 ) if (oi+o 2 ) >D sl 2 /2. 



(31) 

We can now determine 7 self-consistcntly. Given that 
we have an ensemble of TV particles whose rightmost po- 
sition is x max (= In TV), the integrated probability in the 
excess region, J x (Pi(x) + P 2 {x)) dx, should be in the 



order of 1/iV. From Eqs. (23) and (27) we have 

e -ix max = 2v d o 
D sCl N ' 



(32) 



where do is a constant of 0(1). By substituting Eqs. (29) 
and (30) into Eq. (32) and using 02/01 ~ 1/VN, we 
obtain, for small transition rates (oi < D s 7 2 /2), 



= d 1 - 



2oi 
Drf 



+ O (o 2 ) 



(33) 
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where do is another 0(1) constant. 

In the HMF, the transition rate a\ is the inverse of 
the residence time t r of an oscillator near the separatrix 
energy, whose width has been shown to scale as 1/\N. 
We now compute the scaling behavior of this residence 
time, analyzing more closely the dynamics near the en- 
ergy maximum 9 — = n. Consider a particle with 
energy 2-Mjv and phase space coordinates p — and 
9 — 4>n + 7T. By expanding Eq. Q around the potential 
energy maximum, we obtain 



p^ 
2 



2M 



N 



M N 
2 



4>n - T^Y 



(34) 



and the following equations for the single particle dynam- 
ics 



p ~ M N 



5 N 



(35) 



We already know from Sec. |II B| that the global phase c/>n 
exhibits a drift A<fi = 4>N(i) — 4>n(0) — on timescales 
smaller than t<hb ~ N. By integration, we find that 
AO = 9{t) — 0(0) ~ ujt 3 which dominates the dynamics 
of the single particle energy for large times, as 



Ah ~ AO 2 



(36) 



By finally recalling that u> ~ 1/viV, we can determine 
the scaling of the time needed for Aft, to grow up to order 
1/VN, that is 



1 

Oil 



N l/12 i 



(37) 



which is in agreement with numerical results reported in 
Pig. ©a). 

As a result, a± goes to zero algebraically (albeit with 
a small exponent) in the thermodynamic limit and this 
guarantees that the inequality oti < D s j 2 /2 is always 
satisfied asymptotically. By further imposing that x max 
equals to the box width IniV, we have from Eq. (33) 



1 



7 = 



O 



1 



In TV 



(38) 



By now substituting Eq. ( 38 ) into Eq. ( 29 ) , we finally 
obtain 



4 UniV 



(39) 



In other words, because of the divergence of the residence 
time near the separatrix, which induces sustained fluctu- 
ations of the finite-time LE, the asymptotic value of the 
velocity, or the (time-averaged) LE is finite. From the 
quantitative point of view, it is important to notice that 
our estimation for the asymptotic LE 



A, 



A 

4 



(40) 




FIG. 10: (color online), (a) Average single-particle residence 
time t r near the separatrix vs N for U = 0.5. The residence 
time is estimated by measuring the crossing time needed to 
a particle to pass from an energy eL to an energy cr or vice 
versa, where eL < eR are the two energies corresponding to 
the half maximum of the curve Xo(h) reported in Fig. [7] The 
average is taken over 5 x 10 5 to 2 x 10 6 events for each size 
N. The dashed red line indicates N 1 ^ 12 . (b) Single-particle 
finite-time diffusion D(t) vs. time t for U = 0.5. Dashed 
(black), dot-dashed (red) and solid (blue) lines correspond 
corresponds to N = 10 3 , iV = 10 4 , and iV = 10 5 , respectively. 
Inset: Effective diffusion coefficient D s (see text) as a function 
of 1/lniV. Black circles and red squares correspond to U = 
0.5 and U = 0.7, respectively. The dashed lines mark the 
linear extrapolation to the asymptotic value. 



should be interpreted as a lower bound for the actual LE. 
In fact, it has been obtained by assuming the vanishing 
mean velocity for both of the populations, as well as the 
vanishing fluctuations for the still one. Accordingly, there 
are reasons to expect that the contribution to the LE 
from the coupling term may be larger than D s /4. 

In order to compare our estimates with the asymptotic 
values Aoo = 0.056(6) (U = 0.5) and = 0.046(3) 
(U = 0.7) (see Sec. Ill), we need to estimate the effective 



diffusion coefficient D s for a single forced oscillator lying 
at the separatrix energy h = e s . 

From time series ~ 10 6 time units long, we deter- 
mine the mean-square displacement of the integrated 
Lyapunov exponent and thereby, after dividing by the 
elapsed time t, the finite-time diffusion coefficient D{t) 
that is shown in Fig. [To^b) for U = 0.5 and three 
different numbers N of oscillators (which contribute to 
the magnetization). There we notice that, upon increas- 
ing N, D(t) exhibits increasing oscillations of increasing 
period. This is a manifestation of the presence of long 
stretches of positive (negative) local exponents when the 
oscillator is located close to the saddle. In fact, the pe- 
riod is proportional to t s w In N. On the other hand, the 
effective diffusion coefficient D s should be estimated on a 
time scale of the order of the residence time t r close to the 
separatrix and, since we have just seen that t r rj A 1 / 12 , 
it turns out that in the thermodynamic limit t r » t s . 
In the lower-bound spirit of our estimates, we choose to 
identify D s with the minimum of the finite-time diffu- 
sion D(t) in the time interval t € [0, 5t r ]. The results 
are shown in the inset of Fig. [l0|b) for both U = 0.5 
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FIG. 11: (color online). Critical scaling of the largest LE 
Ai (U, N) . (a) Ai as a function of 1/ In N for different energies 
U = 0.70,0.73, and 0.745 from top to bottom, (b) Extrapo- 
lated values of the asymptotic largest LE A^ as a function of 
the distance from the criticality \U— U c \. (c) Size dependence 
of Ai at the critical point U = U c - 



and U = 0.7. By varying the number N of forcing oscil- 
lators from 10 3 to 10 5 and assuming again 1/lniV cor- 
rections, we obtain the asymptotic estimates D s = 0.12 
for U = 0.5 and D s = 0.08 for U = 0.7. It turns out 
that there is approximately a factor two between D s /4 
and Aoo. The main interest of the formula (39) is, how- 



ever, that, representing a lower bound, it shows that the 
largest LE remains strictly positive in the infinite-size 
limit. 



gest that the largest LE exhibits the same critical scaling 
as the magnetization with respect to the system energy 
U — U c , namely, 



Aoo(£0~ \U-U C \^ 2 , for U<U C 



(41) 



In this context Firpo's Riemannian theory [5] predicted 
a different power law A oc (?7) ~ \U - t/ c | 1/6 for U < U c , 
though it was derived under assumptions that are not 
valid near the critical point. 

At criticality, the logarithmic dependence of A on N 
is replaced by an algebraic decay, Ai(f/ c ,iV) ~ N~ x / & , 
toward a vanishing [Fig.|ll[c)]. In fact, this behavior 
can also be explained by the random matrix argument for 
the power-law decay Ai ~ iV -1 / 3 in the disordered phase 
[21 [6]. In the latter case, the disorder of the matrices is 
due to the statistical fluctuations of the magnetization, 
and thus its amplitude 77 scales as 1/vN and the largest 
LE Ai ~ ?7 _2 / 3 ~ TV -1 / 3 . By contrast, at the critical 
point, the disorder is due to the critical decay of the 
magnetization M(U C ,N) ~ N~^l v (see Sec.|lIB), which 
yields 



X(U C ,N) ~ N~ 2 ^ 3 ". 



(42) 



By recalling j3 = 1/2 and v = 2, this indicates A(J7 c , N) ~ 
N^ 1 ^, in agreement with the numerical observation in 
Fig.^c). 



VI. CRITICAL BEHAVIOR OF THE LARGEST 
LYAPUNOV EXPONENT 



So far we have shown that, in the ordered phase 
U < U c , the coupling pressure due to oscillators near 
the separatrix keeps the largest LE Ai positive even in 
the infinite-size limit. This argument does not hold in 
the disordered phase U > U c , where the magnetization 
M (and the separatrix) vanish. Instead, as already men- 
tioned, the largest LE decays to zero as A ~ iV -1 / 3 . 

An interesting question arises then quite naturally: 
what is the behavior of the largest LE in the vicinity of 
the critical point U c l The critical behavior of the largest 
LE A may provide a connection between a dynamic quan- 
tity of the full system and macroscopic thermodynamic 
properties. 

In this section we numerically investigate the criti- 
cal properties of the largest LE Xi(U,N) in the HMF 
model, providing a theoretical account for the observed 
finite-size scaling. We have already seen that, in the or- 
dered phase, Ai decreases logarithmically with increasing 
system size, toward a strictly positive asymptotic value 
Xoo(U). While approaching the critical point, however, 
we find that the logarithmic decay sets in at larger and 
larger sizes [Fig. [Tjja)] and converges to smaller values 
of Aoo (£7). Although large finite-size effects as well as 
critical slowing down prevent us from estimating Aoo(t/) 
near the critical point, our estimates in Fig. |ll[ b) sug- 



VII. THE FULL LYAPUNOV SPECTRUM 

In this section we study the Lyapunov spectrum of the 
HMF model in the ordered phase. This analysis pro- 
vides a more detailed characterization of the instability. 
In particular it helps to assess the (non)extensivity of 
the chaotic dynamics. Given the difficulty of extending 
the theoretical arguments in Sec. [V] beyond the largest 
LE, we restrict our studies to a careful numerical anal- 
ysis. Given the symmetry of the Lyapunov spectrum in 
Hamiltonian systems, it is sufficient to compute the first 
half. 

In Ref. [9] it has been argued that, for U = 0.1, the 
full spectrum vanishes roughly as iV -1 / 3 . However, the 
intermittent behavior observed at low energies [see, e.g., 
Fig. Blb-d)] indicates that the burst state should eventu- 
ally (for N large enough) dominate and thus the TV -1 / 3 
law eventually breaks down. Given this difficulty of deal- 
ing with low energy values, we focus here on a larger- 
energy value, namely U = 0.7. 

Figure |12[ a) shows the Lyapunov spectra Ai as func- 
tions of the rescaled index r = (i — 0.5) /N for different 
system sizes N. This suggests that the spectrum is com- 
posed of two parts: the bulk of the spectrum which de- 
cays toward zero for increasing N, and the initial part 
pinned close to the largest LE, which is clearly visible in 
the inset of Fig. 12 a). 

This scenario is actually coherent with the coexistence 
of extensive and subextensive chaos, recently discovered 
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FIG. 12: Full Lyapunov spectrum at U — 0.7. (a) Lyapunov spectrum Ai as a function of the rescaled index r = (i — 0.5)/N for 
system sizes TV = 32, 64, • • ■ , 1024 from top right to bottom left, as indicated by the arrow. Inset: close-up of the first part, with 
sizes TV = 2048,4096, • • ■ , 16384 added, (b) A, vs N at fixed rescaled indices r = 0.2,0.4,0.6, and 0.8 (from top to bottom). 
The dashed lines indicate Xi ~ 1/VTV. (c) af* N \r'), as defined in Eq. (451, plotted as a function of the logarithmically rescaled 
index r' = (i — 1)/ In TV for TV = 1024, 2048, ■ • ■ , 16384. They show reasonable behavior toward the convergence, implying the 
logarithmic size-dependence (|44|) for these subextensive LEs (see text). 



in generic globally-coupled dissipative systems [IB]. In 
such systems, the Lyapunov spectrum is found to be 
asymptotically flat (a specific realization of extensivity) 
but sandwiched between two vanishing fractions of expo- 
nents located at both ends of the spectrum, where differ- 
ent asymptotic values appear. Finite-size analysis then 
revealed that the bulk of the spectrum scales as 



Ai — A + 



TV 



(43) 



where the asymptotic value Ao corresponds to the LE of 
a single dynamic unit forced by the mean field |16j . 

In Fig. [l2|b) , one can appreciate that the spectrum of 
the HMF decays as predicted by Eq. (|43l with A = 0, 



since a single oscillator has only two variables and thus 
cannot be chaotic. Notice that the existence of these 
zero-Lyapunov bulk components is consistent also with 
the theoretical prediction of Ref . [10] , which did not ex- 
clude the presence of a vanishing (subextensive) fraction 
of different exponents. As for the power-law decay of the 



bulk exponents, the data for small r-values in Fig. 12 'b) 



seem to decrease more slowly than 1/y/N, but this is pre- 
sumably due to strong finite-size corrections induced by 
the bending near the beginning of the spectrum. 

Concerning the subextensive LEs, in dissipative sys- 
tems it was found that there are O(lnTV) exponents 
whose values vary as AW ~ Aoo + a(r') / In TV + C(ln 2 TV) 
with Aoo 7^ A independent of r' , when one fixes a 
logarithmically rescaled index r' = (i — l)/(«o + In TV) 
with a constant io [16 ■ In the HMF, we have shown 
both numerical (Fig. [4]) and theoretical (Sec. [v]) evi- 
dence of this logarithmic dependence for the largest LE, 
Ai ~ Aoo + a(0)/ In TV + b{0) / In 2 TV, with Aoo = 0.046(3) 



for U — 0.7. Now, we assume that the same size- 
dependence holds for subsequent LEs like in dissipative 
systems, with varying coefficients except for the constant 
term: 



Ai ~ Aoo + a(r')/ In TV + b(r')/ In 2 TV. 



(44) 



To examine the validity of this expression, we take the 
Lyapunov spectra \( N \r') = at system size TV and 
compute 



{N) AA (27V) ( r /) m 2 2N _ AX (N) ( r /) ln 2 N 

1 {T)= hi2 ' 



(45) 



with AAW(r') = AW(r') - Aoo- If Eq. M holds, the 
definition in Eq. (45) gives a^ N ' (r') ~ a(r')and the size- 



dependence vanishes. Figure 12 'c) tests this idea and 
indeed verifies that a^(r') approaches an asymptotic 
curve for large sizes TV with the logarithmically rescaled 
index r' = (i— l)/lnTV (here io is set to be zero). There- 
fore, the logarithmic size-dependence (44) holds for these 



subextensive exponents, similarly to dissipative systems. 
Although we need to study larger systems to obtain a 
firmer numerical support, our results on the full spectrum 
of the HMF model are consistent with the coexistence of 
extensive and subextensive exponents, previously found 
for dissipative systems. 



VIII. THE GENERALIZED HMF MODEL 

In order to study the generality of our results, we fi- 
nally turn our attention to a two-dimensional variant of 
the HMF model, introduced by Antoni and Torcini [T5] 
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and later generalized [31 [T7] to the present form. It is denned by the Hamiltonian 
I 



N v 2 +» 2 



1 N 

— {I 1 - cos(xi - Xj)] + [1 - cos(y* - + A [1 - cos(x l - ^) cos(y 4 - %■)]} , (46) 



with two-dimensional coordinates (xi,yi), their conjugate momenta (p x .i,Py,i)y an d a coupling constant A. The 
equations of motion can be written as, 



2/z = J*y,i i 



,4 



2 

,4 



Px,i = -M x sin(x t - <t> x ) - — [P+ sm(xi + y t - ip + ) + P_ sin(xj - y l - ip-)] 
P v ,i = —M y sin(yi - <f> y ) 
with four mean field terms defined as 



[P+ sin(xi + yi - tp + ) - P_ sin(xi - yi - V>-)] > 



i=l t=l 



(47) 



(48) 



As a matter of fact, because of the symmetries of the model, on average M x ~ M y ~ M and P + ~ P_ ~ P and the 
model can be described in terms of two order parameters only. The single oscillator energy can then be written as 

hi = x ' % - v ' 1 + 2 + A - M[cos(x.i - cj> x ) + cos(yi - <p y )] ^-[cos(xi + y t - -0+) + cos(a; l - y l - if>_)] (49) 

(note that in this Section we have not fixed the ground state energy at zero). 



This generalized HMF model is known for its rich and 
generic behavior within the class of systems with long- 
range interactions [3j [17l [18] . For A = it reduces 
to the standard HMF model [Eq. pj]. More generally, 
while the standard HMF model shows a continuous tran- 
sition from the homogeneous to the ferromagnetic, single- 
cluster phase, the generalized HMF model can exhibit 
both continuous and discontinuous canonical transitions 
depending on the value of A, as shown in its phase dia- 
gram (Fig. 13). Moreover, there exists another ordered 
phase, called hereafter the double-cluster phase, which 
is composed of two clusters of oscillators separated on 
average by tt both in Xi and yt, and thus characterized 
by finite values of P and a vanishing magnetization M. 
On the microscopic side, the essential difference with the 
standard HMF model is that here a single oscillator has 
four variables, and hence can be chaotic in the absence of 
any coupling with either an external field or other oscil- 
lators. The largest LE of the full system is therefore not 
purely determined by the coupling effect, unlike in the 
standard HMF model, but receives also a contribution 
from the local dynamics, which depends on the single- 
oscillator energy. 



(point Pi in Fig. 131, in the single-cluster phase, close 



Figure 14 a) shows the largest LE Ai measured for 
three different sets of parameter values: A = 0.2, U = 1.4 



to a (canonical) continuous transition [black circles in 
Fig. ^a,)]; A = 1.0, U = 1.5 (point P 2 ), in the single- 
cluster phase, close to a (canonical) discontinuous tran- 
sition (green diamonds); and A = 6.0, U — 5.0 (point 
P 3 ), in the double-cluster phase (red squares). In all 
the three cases, the largest LE shows the 1/lniV scaling 
toward nonzero asymptotic values, similarly to the stan- 
dard HMF model. From the figure, one notices that in 
proximity of continuous transitions Ai decreases with TV, 
while at Pi , in proximity of the discontinuous transition, 
the maximal LE increases with the system size. This 
peculiarity deserves further investigations. 

Analogously to the ID case, it is instructive to start 
comparing with the behavior of a single oscillator forced 
by constant order parameters. At variance with the ID 
case, the resulting value of the LE can be positive here 
and depends on the energy. When one compares the 
extrapolated asymptotic values Aoo of the full-system LE 
with the maximum value of the single-oscillator LE, Am, 
over possible energy values, we obtain Aoo = 0.16 and 
A M = 0.11 at Pi, Aoo = 0.38 and X M = 0.27 at P 2 , and 
Aoo = 0.23 and Am = at P3; the first asymptotic LE of 
the full system is systematically larger than Am- This is 
again a manifestation of the coupling pressure discussed 
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FIG. 13: (color online). Phase diagram of the generalized 
HMF model. The shaded area denotes the region where mi- 
crocanonical and canonical ensembles differ from one another 
(i.e. the coexistence region of two different thermodynamic 
phases) occurring in correspondence of discontinuous canoni- 
cal transitions. Solid and dashed lines correspond to contin- 
uous and discontinuous transitions, respectively, within the 
canonical ensemble. The red dots indicate the three parame- 
ter values we have studied. 



in Sec. [V] For the first two cases, because of the chaotic 
dynamics of the single oscillator, one does not need to 
introduce the two-family approach taken for the standard 
HMF model, but it is more meaningful to refer to the 
treatment developed in Ref. [16 for dissipative systems, 
which predicts, 



Aoo — + 



D 



(50) 



where D is the diffusion coefficient for the fluctuations of 
the single-oscillator finite-time LE and Ao is the single- 
oscillator LE with an appropriate averaging over energy 
values. By taking into account this correction and us- 
ing A m instead of Ao for the sake of simplicity, we find: 
Aoo « 0.18 (D = 0.15) at P x , and A^ w 0.33 (D = 0.12) 
at Pi- The new values are much closer to the extrapo- 
lated ones, though there is still a remaining gap (espe- 
cially in the second case) , which is presumably due to the 
fact that Eq. (50) was derived under the assumption of 



short ranged time correlations. The slow diffusion across 
different energy surfaces makes this assumption at least 
questionable. In contrast to these cases for the single- 
cluster phase, the situation in the double-cluster phase is 
rather analogous to the standard HMF model; because 
M = in the infinite-size limit, the equations of motion 
in this limit reduce to 



-APsin(ij ±Ui—ip±), 



(51) 



which are equivalent to two uncoupled standard HMF 
models [Eq. pj)]. However, our theoretical approch in 



Sec. [V] should not be applied directly to this case, be- 
cause it deals with finite sizes, where the two variables in 



Eq. (51) are coupled in a non-trivial manner. 



We also studied the full Lyapunov spectrum Aj of the 
generalized HMF model. The results shown in Fig.[l4|b) 
are obtained at Pz- They indicate that the full spectrum 
becomes flatter and flatter for larger sizes, with an ap- 
parent power-law decay of A^ with fixed rescaled index 
r (inset). While the convergence toward zero was ex- 
pected in the ID case, this behavior is questionable for 
the 2D model since in this case the single oscillator may 
be chaotic in the presence of a constant magnetization. In 
fact, it is reasonable to expect that, analogously to the 
dissipative mean-field models discussed in Ref [16], the 
exponents in the bulk of the spectrum converge to the 
value of the LE of a single forced oscillator without cou- 
pling in tangent space. However in a Hamiltonian model 
such as the 2D HMF, it is not clear which exponent one 
should refer to, as it depends on the energy and, more- 
over, the phase space is filled with stable islands. Direct 
measurement of the energy of a single oscillator in the 
full system of size N indicates that, within sufficiently 
long time scales, the single-oscillator energy diffuses as 
largely irrespective of N [Fig. [l4|c)]. Given this exis- 
tence of a well-defined distribution function p(h) for the 
single-oscillator energy, an appropriate reference value Ao 
for the single-oscillator LE would be simply the exponent 
averaged with this distribution function, namely, 



Ao = / dhp{h)X {h) 



(52) 



where \o(h) is the energy-dependent single-oscillator LE. 
The data reported 
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indi- 



m the inset of Fig. 
cate that in the thermodynamic limit Xo(h) vanishes for 
h < es = 2 + A(P + 1) (eg being the single oscillator 
saddle energy in the mean field limit), while finite contri- 
butions arise for larger energy values. At finite N, Ao is 
nothing but the conventional time-averaged LE of a sin- 
gle forced oscillator and is reported in Fig. 14 [d). This 



substantially decreases with increasing N and, in partic- 
ular, in the infinite-size limit, it can reach a positive but 
quite small value of the order of 10 -3 [estimated by a lin- 
ear fit in Fig. |14[d)]. This indicates that the decreasing 
bulk exponents reported in the inset of Fig. [l4|b) can 
have such a small but positive asymptotic value, which 
is however indistinguishable from zero from the available 
numerical data. Moreover, one should notice that each 
oscillator has two nonnegative LEs; the first one can be 
positive or zero as already discussed, while the second one 
is always zero because of the continuous time. It implies 
that one may even expect the occurrence of two bands 
in the asymptotic bulk spectrum, in correspondence with 
these two single-oscillator LEs. Further studies are nec- 
essary to clarify these issues, to elucidate the generality 
of the extensivity and subextensivity found for the stan- 
dard HMF model. 
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FIG. 14: (color online). Lyapunov exponents in the general- 
ized HMF model, (a) Largest LE of the full system, Ai vs 
1/miV for A = 0.2,U= 1.4 (black circles); A = 1.0, U = 1.5 
(green diamonds); and A = 6.0,(7 = 5.0 (red squares). The 
dashed lines indicate linear fits to the scaling regime, (b) Full 
spectrum Ai vs r = (i — 0.5) /N with N = 32,64,128,256, 
and 512 (from top to bottom) for A — 1.0, U = 1.5. Inset: 
Ai vs N for fixed rescaled indices r = 0.5, 1.0, 1.5 from top to 
bottom, (c) Distribution of the single-oscillator energy h for 
A = 1.0, U = 1.5 and systems of size N = 40, 100,400, 1000, 
10000. The arrow indicates increasing system size. In the 
inset: energy-dependent single oscillator LE for the same os- 
cillator numbers N. The vertical dashed line marks the po- 
sition of the mean field saddle energy es = 2 + A(P + 1), 
with the numerical estimate P = 0.469(1). (d) The LE Ao 
of a single oscillator forced by the full system of size N for 
A = 1.0, U — 1.5. The red dashed and green solid lines show 
the results of a quadratic and a linear fitting, respectively, 
which result in finite asymptotic values of Aq. 



IX. CONCLUSIONS AND OPEN PROBLEMS 

The question whether the largest LE in the HMF 
model remains positive or converges to zero in the ther- 
modynamic limit has remained unsettled for a long time. 
We have shown here that the largest LE is indeed positive 
by making use of several subtle properties of globally- 
coupled systems. The first ingredient is what we call the 
"coupling pressure" which induces a finite increase in the 
largest LE (with respect to the LE of a single oscillator 
forced by the mean field). Coupling pressure is a general 
phenomenon that occurs in globally-coupled models of 
both dissipative and conservative dynamical systems and 
arises from the fluctuations of single-oscillator finite-time 
LEs. However, the ID HMF dynamics is even more sub- 
tle, since the LE of the single oscillator under a constant 
field is strictly zero, and its relevant fluctuations must 
be computed by referring to a special type of trajecto- 



ries that come close to the separatrix. More "natural" 
is the behavior of the 2D HMF, since the single oscilla- 
tor dynamics is chaotic and thus the overall scenario is 
analogous to that of standard dissipative chaotic systems 
(see Ref. [TS]). Altogether, our results indicate that the 
thermodynamic limit is rather singular. If one first takes 
the limit N — ¥ oo, no fluctuations can be expected and 
no signature of chaos detected. On the other hand, we 
have shown that the largest LE of an arbitrarily large 
system is always positive. This means that representa- 
tions of the dynamics such as that built in the Vlasov 
equation (which corresponds to assuming N — oo) lose 
the chaoticity of the original dynamics captured by the 
largest LE. 

On a more quantitative level, we have been able to 
derive an explicit expression for a lower bound of the 
largest LE. It would be interesting to improve the argu- 
ment to determine a more accurate estimate and possibly 
predict the dependence on the energy (or, equivalently, 
the temperature). Our numerical analysis indicates that 
the largest LE stays indeed positive in the ordered phase. 
Possibly, it can remain positive in the limit U — > 0, but a 
purely numerical approach is out of question because one 
would need to simulate large enough systems to guaran- 
tee the presence of some oscillators near the saddle of the 
corresponding potential. Near U = 0, the probability for 
an oscillator to come close to the separatrix goes to zero, 
and thus simulations are utterly unfeasible. 

Although we have not been able to extend the theo- 
retical arguments beyond the largest exponent, we have 
undertaken also a general investigation of the entire Lya- 
punov spectrum to investigate the extensivity of the 
chaotic dynamics. Our numerical analysis suggests that 
the asymptotic number of unstable directions is not ex- 
tensive (it grows probably like IniV). It would be de- 
sirable to develop some even approximate argument to 
justify this scaling behavior, which is, so far, only based 
on numerical simulations. 

Finally, we have also analyzed the Lyapunov spectrum 
for the 2D HMF. Such a system is less pathologic than the 
ID model, since the single oscillator dynamics is chaotic 
and it is therefore obvious to expect positive LEs. How- 
ever, a problem remains to be settled regarding the scal- 
ing behavior of the full spectrum. On the basis of all 
arguments developed here and in Ref. [16] , we would ex- 
pect that the bulk of the Lyapunov spectrum (at least 
for r < 1) converges to a finite value. However this is 
not yet seen in our simulations. We cannot exclude that 
this is because the finite value associated to the single 
oscillator dynamics is really small. 

All in all, the results presented here need to be put 
of firmer ground by more rigorous mathematical ap- 
proaches, especially since we have shown that strong 
finite-size effects are at play. The open questions men- 
tioned above also require further work. It is our hope 
that the rather subtle phenomena uncovered here will 
attract such needed attention in the future. 
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